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ABSTRACT 

In  spite  of  extensive  prior  research  on  generation  of  seismic  waves  by  underground  nuclear  explosions,  it  is  still  not 
possible  to  provide  a  complete  explanation  for  the  observed  wavefields,  particularly  at  regional  distances. 

Spherically  symmetric  explosion  models  embedded  in  layered  elastic  media  effectively  model  the  P  phases 
generated  by  explosions,  and  the  major  characteristics  of  some  reflected  and  transmitted  phases.  Nonlinear 
axisymmetrie  Finite  difference  calculations  of  explosions  ineluding  gravity  and  the  effect  of  the  free  surface  can 
model  a  more  realistic  explosion  source  that  directly  generates  shear  w  aves.  These  models  explain  more 
eharaeteristies  of  explosion-generated  seismic  waves,  including  some  aspects  of  regional  shear  phases.  However,  it 
is  clear  that  linear  and  nonlinear  near-souree  3D  effects  are  important  in  many  eases.  SI  I  waves  arc  commonly 
observed  within  a  fewr  km  of  explosions,  too  close  to  have  been  generated  by  (simple)  conversion  of  vertical  and 
radial  components,  and  often  larger  than  those  components.  Furthermore,  it  has  not  been  established  what  impact  3D 
effects  have  on  discriminants  and  on  explosion  yield  estimates.  It  is  important,  therefore,  to  be  able  to  model  and 
understand  how  3D  source  and  source  region  heterogeneity  affect  the  seismic  wavefield,  and  what  impact  this  has 
on  parameters  used  for  nuclear  monitoring. 

We  are  in  the  second  year  of  a  project  to  develop  and  test  a  three-dimensional  nonlinear  Finite  element  code 
CRAM3D,  which  will  be  used  to  calculate  nonlinear  explosion  sources  that  have  both  3D  source  geometry  and  may 
occur  in  a  3D  heterogeneous  medium.  The  code  includes  the  same  well-tested  material  models  that  have  been  used 
in  earlier  axisymmetrie  calculations.  In  addition,  we  are  developing  algorithms  based  on  the  representation  theorem 
to  propagate  the  motion  from  these  source  region  calculations  to  any  desired  distance.  We  have  implemented  a 
technique  that  allows  us  to  propagate  the  results  of  near  source  3D  Finite  element  calculations  to  regional  and 
teleseismie  distances.  The  Green's  function  and  its  derivatives  are  used  in  conjunction  with  the  numerical  solutions 
on  a  monitoring  surface  enclosing  the  complex  source  region.  Full-waveform  solutions  at  distance,  due  to  complex 
explosion  sources,  are  computed  with  the  full- waveform  Green’s  function  using  wavenumber  integration;  surface 
wave  solutions  are  computed  with  the  surface  wave  Green’s  functions  using  mode  summation;  and  far  Field  body 
wave  solutions  are  computed  with  the  outgoing  waves  from  the  source  region.  The  excellent  agreement  in  the 
surface  w  ave  portion  between  the  full-wave  solutions  and  surfaee-wave  solutions  demonstrates  the  accuracy  of  the 
implementation  of  the  representation  theorem  and  the  respective  Green’s  functions  and  their  derivatives. 

To  test  the  code,  we  have  performed  calculations  using  cavities  of  three  shapes:  spherical,  rectangular  and  elliptical, 
each  with  the  same  volume.  An  explosion  with  the  same  yield  was  detonated  inside  each  cavity.  We  compare  the 
solutions  from  these  three  cavity  explosions  in  the  near  Field  and  at  distance.  Gravity  is  included  in  the  calculations, 
and  we  start  with  an  equilibrium  solution  obtained  by  running  the  finite-element  CRAM3D  with  overburden 
pressure  only,  prior  to  the  start  of  the  explosion  calculation.  Nonlinear  deformation  is  seen  around  the  cavity.  The 
results  show  very  good  agreement  between  2D  and  3D  solutions  at  distance  for  the  spherical  cavity  explosion. 
Nonspherieal  wave  components  from  nonspherieal  rectangular  and  elliptical  cavities  are  clearly  seen  in  the  near 
Field  The  rectangular  cavity  show  s  more  pronounced  tangential  motion  than  the  elliptical  cavity  away  from  axes  of 
symmetry. 
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OBJECTIVES 

The  objective  of  this  project  is  to  investigate  the  generation  of  complex  seismic  waves  by  explosions  in  media  with 
3D  heterogeneity  using  a  method  based  on  the  exact  representation  theorem  for  propagating  complex  3D  source 
calculations  to  local,  regional  and  teleseismic  distances. 

RESEARCH  ACCOMPLISHED 


Introduction 

Source  physics,  near  source  scattering,  and  propagation  effects  are  all  important  to  understanding  seismic  phases 
used  in  nuclear  monitoring.  Significant  bodies  of  literature  exist  that  address  each  subject  individually.  In  addition  to 
extending  our  understanding  of  source  physics  and  near  source  scattering  to  include  the  effects  of  realistic  3D 
heterogeneity,  this  project  links  the  progression  of  energy  from  its  generation  by  the  source,  through  the  near  source 
region,  and  into  its  partitioning  among  local,  regional,  and  teleseismic  phases.  Distinguishing  the  far-ficld  P-  and 
S-wavcs  enables  us  to  quantify  the  effects  of  3D  structure  on  both  P-  and  S-wave  generation.  Complete  regional 
waveforms  show  how  this  energy  is  partitioned  among  the  distinct  phases,  which  are  important  to  event  detection, 
identification,  and  magnitude  estimation.  How  energy  is  distributed  among  surface  wave  modes  determines  Lg 
amplitudes  and  Rg  amplitudes,  including  near  source  conversions  between  these  phases.  Modal  excitation  of  Lg  as 
well  as  Rg  has  a  significant  depth  dependence  that  is  often  neglected  in  nuclear  monitoring  studies  (Baker  ct  al., 
2004),  and  these  amplitudes  can  be  disrupted  further  by  near  source  effects  and  source  region  structure.  This  can 
either  degrade  or  improve  discrimination  capability  depending  on  how  well  it  is  understood. 

While  it  is  not  surprising  to  observe  SH  waves  from  any  single  event,  it  is  surprising  to  note  it  is  always  present;  at 
regional  distances,  SII  scales  with  yield  about  as  well  as  P-waves  While  there  arc  many  mechanisms,  such  as 
ncar-source  scattering,  tectonic  release,  etc.  that  can  generate  Si!  waves,  all  of  these  effects  should  be  highly 
variable  from  one  event  to  the  next.  To  address  this  question,  wc  need  to  look  at  the  types  of  3D  source  effects  that 
can  exist,  the  range  of  variability  that  would  be  expected  from  them,  and  w  hether  the  predictions  are  consistent  with 
observations. 

Souree  Region  Calculations  and  Propagation  Using  the  Representation  Theorem 

Our  approach  is  to  perform  3D  explosion  source  region  calculations,  and  then  to  propagate  the  wavcficld  to  local, 
regional  and  teleseismic  distances  using  layered  earth  Green’s  functions.  Wc  are  interested  in  near-source 
heterogeneities  in  both  the  nonlinear  and  linear  regimes,  and  therefore  require  both  nonlinear  and  linear  3D  codes  to 
model  the  source  region.  In  previous  projects,  we  have  used  two  nonlinear  codes,  STELLAR  and  CRAM,  which  are 
described  briefly  in  Table  1  (we  have  also  used  the  1 D  nonlinear  code  SKIPPER,  which  is  a  spherically  symmetric 
version  of  CRAM),  and  the  2D  and  3D  linear  elastic  code  TRES3D. 


Table  1.  Numerical  simulation  tools  used  in  this  project 


Numerical  Simulation  Tools 

STELLAR 

Eulerian  finite  difference  code.  Used  to  simulate  the  early  time 
history  of  the  explosion  shock.  It  handles  material  strength  correctly, 
which  is  difficult  for  a  Eulerian  code.  Uses  second  order  accurate 
Riemann  solver  scheme.  1 D,  2D  planar  and  axisymmetnc,  and  3D. 

CRAM 

Lagrangian  nonlinear  finite  difference  code.  Has  been  used 
extensively  for  axisymmetric  explosion  calculations.  A  3D  version 
of  the  code  is  being  developed  in  this  project. 

TRES3D 

Elastic  finite  difference  code.  2D  planar  and  axisymmetric,  and  3D. 

Elastodynamie  Representation  Theorem 

The  time  dependent  displacements  and  stresses  from  3D  source 
region  calculations  are  saved  on  a  monitoring  surface  localed  outside 
of  the  region  of  nonlinear  response  and/or  3D  heterogeneity.  A 
numerical  implementation  of  the  representation  integral  is  then  used 
to  compute  the  corresponding  far-field  seismic  radiation 

In  past  projects,  we  have  used  these  codes  in  the  following  ways: 
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1 .  Axisymmetrie  CRAM  was  used  together  with  the  representation  theorem  to  propagate  the  results  of 
nonlinear  axisymmetrie  finite  difference  calculations  to  regional  and  teleseismie  distances  (Day  et  al., 

1987;  Stevens  et  al.,  1991;  Stevens  et  al,  2004). 

2.  Axisymmetrie  STELLAR  was  used  to  calculate  the  early  stages  of  an  explosion.  The  solution  at  an 
appropriate  time  was  overlain  onto  CRAM  to  be  propagated  out  to  the  linear,  elastic  region.  The 
representation  theorem  was  used  to  propagate  the  waveform  to  regional  and  teleseismie  distances  (Rimer  et 
al,  1994).  A  similar  technique  was  used  to  overlay  results  from  the  Los  Alamos  National  Laboratory  code 
SOIL  onto  CRAM  and  propagate  the  results  (Davis  et  al,  1992). 

3.  3D  STELLAR  was  used  to  calculate  near-field  waveforms  from  explosions  in  rectangular  cavities  (Stevens 
et  al,  2006). 

4  TRES3D  was  used  to  calculate  the  scattering  from  explosions  in  a  region  with  3D  heterogeneity  and 
topography  (Stevens  et  al,  2004). 

In  the  current  project,  we  are  doing  the  following:  1 )  using  STELLAR  to  perform  very  near  source  nonlinear  3D 
calculations  and  TRES3D  to  perform  linear  clastic  near  source  and  source  region  3D  calculations;  2)  developing  a 
3D  version  of  CRAM  to  model  explosions  from  the  source  out  through  the  very  important  nonlinear  to  linear 
transition  region;  and  3)  completing  the  implementation  of  the  elastodynamic  representation  theorem  for  full 
waveforms,  modes  and  body  waves. 

The  technique  for  propagating  numerical  calculations  using  the  representation  theorem  is  to  save  displacements  and 
stresses  on  a  monitoring  surface  surrounding  the  nonlinear  and/or  heterogeneous  region  of  the  calculation,  and  then 
to  convolve  these  with  a  Green’s  function  appropriate  for  the  external  region  (Stevens  et  al,  1991).  In  the  eases  that 
we  have  done  previously,  two-dimensional  axisymmetrie  nonlinear  finite  difference  calculations  were  performed  to 
model  the  nuclear  explosion,  and  the  stresses  and  displacements  from  the  calculation  were  saved  on  a  cylindrical 
surface  in  the  elastic  region  outside  the  region  of  complex  nonlinear  behavior.  We  then  invoked  the  representation 
theorem  and  integrated  the  stresses  and  displacements  with  an  axisymmetrie  Green’s  function  to  calculate  the 
displacement  at  any  point  outside  of  the  calculation.  We  performed  such  calculations  in  2D,  using  Green’s  functions 
for  far-field  body  waves,  for  modes  and  for  full  regional  waveforms  using  wavenumber  integration.  The  equations 
for  the  Green’s  functions  for  surface  waves  are  given  by  Baehe  et  al.  (1982).  The  Green’s  functions  for  the  complete 
seismograms  are  computed  using  a  ring  load  source,  from  an  algorithm  based  on  the  work  of  Lueo  and  Apsel  (1983) 
and  Apsel  and  Lueo  (1983).  The  Green’s  functions  for  body  waves  are  generated  by  a  procedure  similar  to  that 
described  by  Baehe  and  Harkrider  (1976)  using  a  saddle  point  approximation  to  calculate  a  far-field  plane  wave  for 
a  given  takeoff  angle  from  a  source  in  a  plane-layered  medium.  Our  objective  in  using  multiple  types  of  Green’s 
functions  is  to  gain  as  much  insight  as  possible  into  the  nature  of  the  seismic  wavefield  generated  by  the  source  An 
important  part  of  the  current  project  is  to  adapt  these  techniques  to  propagate  seismic  waves  from  3D  source 
calculations.  Although  any  closed  surface  can  be  used  for  representation  theorem  integration,  we  use  a  cylindrical 
surface  for  axisymmetrie  problems,  and  a  rectangular  surface  for  3D  problems. 

3D  implementation  of  the  representation  theorem 

The  key  to  extending  the  axisymmetrie  representation  theorem  discussed  above  to  3D  is  to  recognize  that  while  the 
deformation  in  the  source  region  may  be  arbitrarily  complex,  if  the  structure  can  be  approximated  as  a  plane-layered 
medium  outside  of  the  source  region,  then  the  knowu  Green’s  function  for  a  plane-layered  medium  applies  (this  also 
assumes  that  we  can  neglect  the  interaction  of  any  baekseattered  waves  returning  to  the  source  region  after  leaving 
it).  Note  that  the  representation  theorem  is  exact.  That  is,  no  matter  how  complex  the  3D  motion  is  on  the  source 
region  boundary,  it  will  be  correctly  propagated  by  the  representation  theorem.  The  following  benchmark  test 
demonstrates  the  performance  of  the  method  by  comparing  the  results  with  equivalent  finite  difference  calculations 
and  wavenumber  integration  seismograms. 

The  representation  theorem  states  that*displaeement  at  an  observation  point  is  made  up  of  contributions  due  to  body 
forces  throughout  the  source  volume,  plus  contributions  due  to  the  traction  and  displacement  on  the  source  volume 
surface  (Aki  and  Richards,  1 980)  In  the  three-dimensional  numerical  finite  difference  calculations,  we  save 
displacements  and  stresses  due  to  the  seismic  source  on  a  monitoring  surface  on  the  boundary  of  a  rectangle  (five 
planar  surfaces,  excluding  the  upper  surface),  and  calculate  Green’s  functions  from  each  point  on  the  monitoring 
surface  to  the  receiver  and  thus,  the  synthetic  seismogram  at  the  receiver  point  X  outside  of  the  monitoring  surface  is 

obtained  by  integrating  over  the  monitoring  surface  su  : 
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", = <#s„  (**)  *  T“  (fl- »;  (o*  (*•*>*  c) 

in  the  frequency  domain,  where  G'  (£;  A")  and  S  (£;  X)  are  the  Green’s  function  and  the  stress  tensor  on  the 

monitoring  surface  due  to  a  unit  impulsive  force  at  X  in  direction  /,  T}  7  is  the  traction  on  the  monitoring  surface  due 

to  the  seismic  source,  u  is  the  displacement  on  the  monitoring  surface,  and  n  is  the  normal  to  the  monitoring  surface. 
The  operator  *  denotes  convolution  and  the  summation  convention  is  assumed 

Equation  1  is  applicable  to  any  Green’s  function  for  the  exterior  model,  and  so  we  can  use  a  full  waveform  Green’s 
function,  far-field  body  wave  Green’s  function,  and/or  modal  Green’s  function  as  we  have  discussed  earlier  for 
axisymmetric  problems.  We  have  implemented  this  technique  for  full  waveform  seismograms  calculating  the 
Green’s  functions  using  wavenumber  integration.  We  have  also  implemented  the  technique  for  body  waves  and  for 
modes  using  corresponding  Green’s  functions. 

The  synthetic  seismograms  are  computed  using  the  following  steps: 

1 .  Displacements  and  stresses  are  saved  on  the  monitoring  surface  during  the  finite  difference  calculations. 

2.  If  necessary,  the  monitonng  solutions  arc  resampled  onto  a  coarse  grid,  as  permitted  by  the  required 
resolution.  If  there  is  a  symmetry  boundary,  the  entire  monitonng  surface  is  constructed  first. 

3.  The  finite  difference  solutions  at  each  point  on  the  monitonng  surface  are  transformed  into  the  frequency 
domain. 

4.  The  displacement  and  stress  Green’s  functions  due  to  the  three  orthogonal  forces  at  the  receiver  location 
are  calculated  for  each  location  on  the  monitoring  surface  in  the  frequency  domain  (using  reciprocity). 

5.  Equation  ( 1 )  is  used  to  obtain  the  solution  at  the  receiver  in  the  frequency  domain. 

6.  The  solution  is  transformed  back  to  the  time  domain. 

Step  5  is  implemented  with  a  generalized  interface  code,  which  takes  as  inputs  the  monitoring  wavefields  and  the 
Green’s  functions  and  stress  fields  for  any  of  the  full-waveforms,  surface  waves  or  body  waves.  Specifically, 
full-wavefonn  solutions  at  distance,  due  to  complex  explosion  sources,  are  computed  with  the  full-waveform 
Green’s  function  using  wavenumber  integration.  Surface  wave  solutions  are  computed  with  only  with  the  surface 
wave  Green’s  functions  using  mode  summation.  Far  field  body  wave  solutions  are  computed  with  the  outgoing 
waves  from  the  source  region.  The  principal  advantage  of  this  approach  is  that  it  allows  us  to  perform  detailed 
calculations  of  the  source  region  and  then  propagate  the  results  to  distances  that  would  be  impractical  or  impossible 
to  include  in  the  same  numerical  calculation.  In  addition  to  reducing  cost  and  time,  the  hybrid  method  is  also  more 
accurate,  as  numerical  dispersion  increases  with  the  size  and  duration  of  numerical  calculations. 

Benchmark  Tests  with  Gravity 

We  have  implemented  an  explicit  three-dimensional  Lagrangian  finite  element  algorithm  that  is  capable  of  using 
multiple  processors  (Stevens  and  Xu,  2008).  All  of  the  nonlinear  material  models  from  2D  CRAM  have  been 
implemented  and  gravity  is  included.  The  cavity  is  placed  near  the  center  of  the  grid  and  is  enclosed  by  a  spider  grid 
which  facilitates  applying  the  pressure  boundary  condition  and  rezoning  elements,  as  implemented  in  the  two- 
dimensional  axisymmetrie  code,  CRAM. 

To  test  the  code,  we  have  performed  calculations  using  three  kinds  of  cavities:  spherical  (radius  5m),  rectangular 
(each  side  8.06m)  and  elliptical  (three  axis  lengths  are  6m,  5m  and  4.1667m),  each  with  the  same  volume.  The  same 
yield  (0.2kt)  explosion  is  detonated  in  each  cavity.  The  material  external  to  the  cavity  is  a  model  for  Degelen  granite 
(Stevens  et  al.,  2003).  The  model  consists  of  two  layers:  the  top  layer  is  30m  thick  elastic  and  the  bottom  layer  is 
nonlinear.  Gravity  is  included  in  the  calculations,  and  we  start  with  an  equilibrium  solution  obtained  by  running  the 
finite-element  CRAM3D  with  overburden  pressure  only,  prior  to  the  start  of  the  explosion  calculation.  We  compare 
the  solutions  from  these  three  cavity  explosions  in  the  near  field  and  at  distance.  The  shapes  of  the  cavities  are 
shown  in  Figure  1  and  the  seismograms  in  the  near  field  are  shown  in  Figure  2.  Each  seismogram  plot  corresponds 
to  a  receiver  in  the  same  pattern  in  Figure  1.  The  radial,  tangential  and  vertical  components  are  represented  by  red, 
green  and  blue  lines.  The  top  right  plots  correspond  to  the  spherical  cavity  and  show  no  visible  tangential  motion 
(green  lines)  at  these  receivers.  The  bottom  left  plots  correspond  to  the  rectangular  cavity  and  show  clear  tangential 
motion  off  the  symmetry  axes  (vertical,  horizontal  and  diagonal).  The  bottom  right  plots  correspond  to  the  elliptical 
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cavity  and  show  some  visible  tangential  motion  off  the  symmetry  axes  (vertical  and  horizontal).  The  main  waveform 
characteristics  in  the  elliptical  case  are  quite  similar  to  those  in  the  spherical  case. 

SPHERE  RECTANGULAR  ELLIPSE 


Figure  1.  Three  cavity  shapes  with  the  same  volume  used  for  3D  nonlinear  explosion  calculations. 
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Figure  2.  Near  field  receiver  locations  (top  left):  blue  circle  indicates  the  cavity  location.  The  near  field 

waveforms  due  to  different  cavities  are  shown  at  top  right  and  bottom.  Red,  green  and  blue  lines 
correspond  to  radial,  tangential  and  vertical  components.  Note  the  tangential  components  from  the 
non-spherieal  cavities  (bottom  two  figures). 
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All  the  three  explosions  yield  nonlinear  deformation  around  the  cavity  in  Figure  3.  The  horizontal  (top)  and  vertical 
(bottom)  slices  across  the  cavity  center  are  shown  for  comparison.  Again,  both  spherical  and  elliptical  cavities  show 
similar  patterns  but  the  nonlinear  deformation  region  is  enlarged  for  the  rectangular  cavity  relative  to  the  other  two 
cavities,  most  likely  caused  by  the  effects  of  the  comers. 
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Figure  3.  Nonlinear  deformation  extent  around  the  cavities.  The  top  plots  are  the  horizontal  slices  across  the 
cavity  center  and  the  bottom  illustrate  vertical  slices.  The  asymmetry  in  the  Z  direction  is  caused  by 
the  variation  of  overburden  pressure  w  ith  depth.  More  nonlinear  deformation  occurs  above  the 
cavity  than  below. 

Far  field  solutions  are  obtained  with  the  representation  theorem  (Equation  1)  by  using  the  displacements  and  stresses 
recorded  on  the  monitoring  surfaces  in  the  finite-element  calculations.  An  interface  code  was  developed  to  correctly 
match  the  numerical  solutions  with  the  Green’s  functions  and  the  corresponding  stresses.  We  compare  the  full 
waveforms  and  surface  waves  at  two  different  locations  with  the  same  distance.  One  is  at  x=2000m,y=  1500m 
(location  1,  no  symmetry)  on  the  surface  and  the  other  at  x=2500m  and  y=0m  (location  2  at  the  symmetry  axis).  The 
full  waveforms  are  computed  by  the  wavenumber  integration  method  and  the  surface  waves  are  computed  using 
mode  summation.  The  results  at  the  two  locations  for  the  three  cavities  are  shown  in  Figure  4.  The  left  corresponds 
to  the  location  1  and  the  right  location  2.  Each  panel  has  three  components.  The  red  and  green  lines  indicate  the  full 
waveform  solutions  and  surface  wave  solutions,  respectively,  and  the  blue  dashed  lines,  the  2D  full  waveform 
solutions  in  the  spherical  cavity  case.  All  the  waveforms  are  low  pass  filtered  below  5Hz.  It  is  clearly  seen  that  the 
mode  summation  solutions  (green  lines,  Figure  4)  match  the  surface  wave  portions  of  the  full  waveform  solutions 
(red  lines,  Figure  4)  for  all  the  cavity  types  at  two  locations  very  well.  For  the  spherical  cavity  full-waveform 
solutions,  there  is  also  very  good  agreement  between  2D  (blue  dashed  lines,  Figure  4)  and  3D,  validating  the  proper 
implementation  of  the  Lagrangian  finite-element  algorithms  and  the  representation  theorem  in  3D.  It  is  also  noted 
that  the  spherical  and  elliptical  cavities  have  the  similar  waveforms  at  distance,  as  in  the  near  field.  The  wave 
amplitudes  arc  slightly  larger  for  the  rectangular  cavity  (center  row.  Figure  4)  and  consistent  with  nonlinear 
deformation  extent.  The  tangential  motions  arc  very  small  for  the  three  cavities  and  indicate  that  the  source 
asymmetry  due  to  a  small  yield,  as  seen  in  the  horizontal  nonlinear  deformation  distribution  (top  row.  Figure  4),  is 
quite  weak  at  low  frequencies,  which  is  also  verified  by  comparing  the  waveform  solutions  at  the  two  locations. 


579 


2009  Monitoring  Research  Review:  Ground-Based  Nuelear  Explosion  Monitoring  Technologies 


Figure  4.  Far  field  waveforms  using  the  wavenumber  integration  method  and  mode  summation.  Waveforms 
at  location  1  are  shown  on  the  left,  location  2  on  the  right.  There  is  excellent  agreement  between  the 
full-waveform  solutions  and  surface  wave  solutions  for  the  surface  wave  part  of  the  waveform. 
Excellent  agreement  is  also  demonstrated  for  the  full  waveform  solutions  in  2D  (blue  dashed  lines) 
and  3D  (top). 
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CONCLUSIONS  AND  RECOMMENDATIONS 

We  are  in  the  seeond  year  of  a  project  to  understand  3D  effects  on  seismic  radiation  from  underground  nuclear 
explosions.  We  have  nearly  completed  development  of  a  3D  version  of  CRAM,  the  Lagrangian  code  we  have  used 
previously  for  performing  axisymmetric  calculations  of  underground  explosions.  We  have  also  implemented  an 
interface  code,  which  utilizes  any  Green’s  functions  in  order  to  propagate  the  results  of  the  near  field  3D 
calculations  to  regional  and  teleseismie  distances  using  the  representation  theorem.  The  objectives  arc  to  complete 
implementation  of  the  numerical  methods  and  then  perform  3D  calculations  to  understand  and  model  the  effects  of 
3D  source  region  heterogeneity  and  the  seismic  response  to  it  in  a  realistic  source  scenario 
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